Programming Assignment 1: Percolation

通过蒙特卡洛模拟方法来估计渗流阈值。

Percolation. 给一个有随机分布的绝缘和金属材料的组成的复合系统。例如我们想知道哪些部分必须是金属材料才能让这个复合系统是一个电导体。或者在一个多孔的地形，在表面有水或者油，在什么情况下水或者油能够从最表面渗透到最底层。科学家把这种过程的模型叫做Percolation。

The model. 在Assignment中，用一个NxN的格子表示percolation系统，每一个格子是打开或者关闭，打开是白色关闭是黑色。如果一个格子是full，首先他必须是打开额，然后表示从最顶上通过相连(4方向)的打开的格子可以渗透到这个位置。当一个系统是percolates，表示能从最顶层渗透到最底层，也就是说，最底层存在打开的格子是full。

 

The Problem. 研究人员对一下的问题感兴趣，如果每一个格子是独立的，并且被打开的概率为p，那么系统percolates的概率是多少？p=0，percolates概率为0，p=100，percolates的概率为100。下图是20x20和100x100格子的概率p的分布：

当N足够大时, 有一个阈值P*, 使得当p < p*时候，任意的N*N网格，几乎不能被渗透, 并且当p > p*, 基本能够被渗透。p*没有准确的数值解。任务是写一个计算估计p*的算法。

Percolation data type.

public class Percolation {
   public Percolation(int N)              // create N-by-N grid, with all sites blocked
   public void open(int i, int j)         // open site (row i, column j) if it is not already
   public boolean isOpen(int i, int j)    // is site (row i, column j) open?
   public boolean isFull(int i, int j)    // is site (row i, column j) full?
   public boolean percolates()            // does the system percolate?
}

其中网格的左上角为(1, 1)，右下角为(N, N)。

使用boolean[] matrix来标记网格的打开和关闭。

isOpen表示当前网格是否是打开，直接返回matrix当前位置的值即可。

思考open，isFull，percolates的使用并查集如何实现和复杂度。

open要做什么？思路比较简单，将一个关闭的网格位置打开，并且和周围4方向上已经打开的格子连通在一起。

isFull判断当前点是否能够渗透到，最简单的想法，枚举第一层每个打开的格子，并查集查询判断与当前格子是否在一个集合中，如果在那么就是成功的。复杂度枚举最差N, 再加上每次需要进行一次1个Find()的判断。

percolates判断系统是否能够被渗透，最简单的想法，枚举第一层和最后一层打开的格子，判断是否在一个集合中，如果在那么渗透成功，复杂度N^2，再加上并查集Find()判断。

最简单的想法太暴力，复杂度太高，那么我们开始优化

1. 加入虚拟节点，首尾加入两个虚拟节点，首虚拟节点与第一层所有打开的元素相关联，尾虚拟节点和最后一层所有打开的节点相关联，那么对于isFull操作，如果首虚拟节点与当前点在一个集合，那么成功，对于percolates操作如果首尾虚拟节点在一个集合中，那么系统是渗透成功的。这样所有的这些操作中都不会有枚举的for循环出现了，但是也会出现一个问题，backwash，因为有一些格子虽然从上面不能渗透到，但和尾虚拟节点连接后，他们也和首虚拟节点在一个集合了。

2. 如何避免backwash？想了个比较方法，又加入了一个并查集，现在两个并查集，一个只负责维护首虚拟节点，当进行isFull判断是否只考虑这个并查集，另一个并查集进行首尾虚拟节点的维护，用在percolates的判断中，这样backwash的问题也解决了，但是牺牲了空间。

下面是实现代码：

public class Percolation {
    
    private boolean[] matrix;
    private int row, col;
    private WeightedQuickUnionUF wquUF;
    private WeightedQuickUnionUF wquUFTop;
    private boolean alreadyPercolates;
    
    public Percolation(int N) {
        if (N < 1) throw new IllegalArgumentException("Illeagal Argument");
        wquUF = new WeightedQuickUnionUF(N*N+2);
        wquUFTop = new WeightedQuickUnionUF(N*N+1);
        alreadyPercolates = false;
        row = N;
        col = N;
        matrix = new boolean[N*N+1];
    }
    
    private void validate(int i, int j) {
        if (i < 1 || i > row) 
            throw new IndexOutOfBoundsException("row index i out of bounds");
        if (j < 1 || j > col) 
            throw new IndexOutOfBoundsException("col index j out of bounds");        
    }
    
    public void open(int i, int j) {
        validate(i, j);
        int curIdx = (i-1)*col + j; 
        matrix[curIdx] = true;
        if (i == 1) {
            wquUF.union(curIdx, 0);
            wquUFTop.union(curIdx, 0);
        }
        if (i == row) {
            wquUF.union(curIdx, row*col+1);
        }
        
        int[] dx = {1, -1, 0, 0};
        int[] dy = {0, 0, 1, -1};
        for (int dir = 0; dir < 4; dir++) {
            int posX = i + dx[dir];
            int posY = j + dy[dir];
            if (posX <= row && posX >= 1 
                    && posY <= row && posY >= 1 
                    && isOpen(posX, posY)) {
                wquUF.union(curIdx, (posX-1)*col+posY);
                wquUFTop.union(curIdx, (posX-1)*col+posY);
            }
        }
    }
    
    public boolean isOpen(int i, int j) {
        validate(i, j);
        return matrix[(i-1)*col + j];
    }
    
    public boolean isFull(int i, int j) {
        validate(i, j);
        int curIdx = (i-1)*col+j;
        if (wquUFTop.find(curIdx) == wquUFTop.find(0)) return true;
        return false;
    }
    
    public boolean percolates() {
        if (alreadyPercolates) return true;
        if (wquUF.find(0) == wquUF.find(row*col+1)) {
            alreadyPercolates = true;
            return true;
        } 
        return false;
    }
    
    public static void main(String[] args) {
        Percolation perc = new Percolation(2);
        perc.open(1, 1);
        perc.open(1, 2);
        perc.open(2, 1);
        System.out.println(perc.percolates());
    }
    
}

Monte Carlo simulation. 估计percolation的阈值，初始化时候格子都是关闭的，随机寻找一个关闭的位置打开，直到系统可以渗透为止，打开的格子比上总格子数就是阈值。

这一部分就是一些数值计算比较简单，不赘述了。

public class PercolationStats {
    private double[] x;
    private int expTimes;
    
    public PercolationStats(int N, int T) {
        if (N <= 0 || T <= 0) 
            throw new IllegalArgumentException("Illeagal Argument");
        x = new double[T+1];
        expTimes = T;
        for (int i = 1; i <= T; ++i) {
            Percolation perc = new Percolation(N);
            boolean[] isEmptySiteLine = new boolean[N+1];
            int numOfLine = 0;
            while (true) {    
                int posX, posY;
                do {
                    posX = StdRandom.uniform(N)+1;
                    posY = StdRandom.uniform(N)+1;
                } while (perc.isOpen(posX, posY));
                perc.open(posX, posY);
                x[i] += 1;
                if (!isEmptySiteLine[posX]) {
                    isEmptySiteLine[posX] = true;
                    numOfLine++;
                }
                if (numOfLine == N) {
                    if (perc.percolates())
                        break;
                }
            }
            x[i] = x[i] / (double) (N*N);
        }
    }
    
    public double mean() {
        double mu = 0.0;
        for (int i = 1; i <= expTimes; ++i) {
            mu += x[i];
        }
        return mu / (double) expTimes;
    }
    
    public double stddev() {
        if (expTimes == 1)
            return Double.NaN;
        double sigma = 0.0;
        double mu = mean();
        for (int i = 1; i <= expTimes; ++i) {
            sigma += (x[i]-mu)*(x[i]-mu);
        }
        return Math.sqrt(sigma / (double) (expTimes-1));
    }
    
    public double confidenceLo() {
        double mu = mean();
        double sigma = stddev();
        return mu - 1.96*sigma / Math.sqrt(expTimes);
    }
    
    public double confidenceHi() {
        double mu = mean();
        double sigma = stddev();
        return mu + 1.96*sigma / Math.sqrt(expTimes);
    }
    
    public static void main(String[] args) {
        int N = Integer.parseInt(args[0]);
        int T = Integer.parseInt(args[1]);
        PercolationStats percStats = new PercolationStats(N, T);
        StdOut.printf("mean = %f\n", percStats.mean());
        StdOut.printf("stddev = %f\n", percStats.stddev());
        StdOut.printf("95%% confidence interval = %f, %f\n", 
                      percStats.confidenceLo(), percStats.confidenceHi());
        
    }
    
    
}